Stability of topological solitons in modified two-component Ginzburg-Landau model 
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We study the stability of Hopfions embedded in a certain modification Ginzburg-Landau model of 
two equally charged condensates. It has been shown by Ward [Phys. Rev. D66, 041701 (R) (2002)] 
that certain modification of the ordinary model results in system which supports stable topological 
solitons (Hopfions) for some values of the parameters of the model. We expand the search for stability 
into previously uninvestigated region of the parameter space, charting an approximate shape for the 
stable/unstable boundary and find that, within the accuracy of the numerical methods used, the 
energy of the stable knot at the boundary is independent of the parameters. 
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I. INTRODUCTION 

Topological solitons, be it vortices, knots, instantons 
or other objects, enjoy widespread interest within many 
fields of physics, perhaps most notably in the fields of par- 
ticle physics and condensed matter, where these objects 
invariably occur as solutions of the field equations. In 
condensed matter physics, topologically stable vortices 
are also a routinely seen in experiments. Therefore, it 
is crucial to understand the basic properties of topolog- 
ical solitons in the Ginzburg-Landau and related mod- 
els. It is against this background that we have studied 
the static Ginzburg-Landau model, also known as the 
Abelian Higgs model. Some years ago, it was demon- 
strated that the archetypal 3D classical field model sup- 
porting topologically stable closed vortices, the Faddeev- 
Skyrme model 0, 0, S, II H H, 0, B H , can be embedded 
in the Ginzburg-Landau model by a change of variables 
. There was also an earlier work, where the FS model 
was obtained from Ginzburg-Landau model in a deriva- 
tive expansion [111 ], but this method does not allow for 
investigation of solitons in Ginzburg-Landau model since 
it does not correspond to any parameter limit. It was 
further conjectured in Ref. |l0j that the two-component 
Ginzburg-Landau model should, due to this embedding, 
also support the same topological structures as the FS 
model does, but more recent studies do not support this 
conjecture [H, [T^]. However, Ward found [ijj that by 
modifying the model suitably, stable closed vortices ap- 
pear as minimum energy configurations of the theory. 
We expand on Ward's work and find the stable/unstable 
boundary in the (k, if) parameter space of the model. 
We also find that the energy of the stable minimum en- 
ergy configuration along the said boundary is constant 
(to within the accuracy of the methods used) . The results 
presented here should also be useful in constructing such 
initial configurations in Ginzburg-Landau model which 
relax into a local, non-zero energy minimum instead of 
the global one and thus provide a way to construct knot- 



ted solitons in the Ginzburg-Landau model. 



II. THE MODEL 

The static Abelian Higgs model with two charged 
Higgs bosons is mathematically the same as the 
Ginzburg-Landau model with two flavors of Cooper pairs 
or super-fluids. The paper will use the following no- 
tations. The indices run as follows: j, k, I £ {1,2,3}, 
a 6 {1,2}, n, v, i £ {0,1,2,3} and the fields are <F = 
(V>i ^) T , = d^A v ~ d v A^ B = e^ kl d k Ai and the 
gauge-covariant derivative has the form D M = d^—ij^A^; 
when working in three dimensions (fi £ {1,2,3}) we 
will also write D = V — *|§^.- With these notations, 
the standard Lagrangian density of the two-component 
Ginzburg-Landau model can be written as 



1m a 



\F^, (1) 



which gives the static energy density 



v(1>wh) + ifa\\B 



(2) 



where we have used SI units. The form of the potential 
is not very important as long as it maintains the SU (2) 
symmetry of >F and enforces the condition = con- 
stant ^ at some limit of the parameters of the potential; 
here we have used 



(3) 



The electric coupling constant displays an explicit factor 
of 2 due to the interpretation of the Ginzburg-Landau 
model as a superconductor, where the %j) a become Cooper 
pairs. For computational purposes, it is practical to use 
natural units, where h = c = /j,q — 1 and rescale the 
fields by tp a — * ^ a ^m a . For flexibility, we retain a freely 
selectable electric charge but replace 2e — * g, finally ob- 
taining the energy density (now = dk — igAk) 



£ = h\Dip 



B\\ 



(4) 
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The remainder of this paper will be in natural units. 
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The embedding of Babaev et al. [10( is such that a 
closed vortex can be defined by the fields ip a , leaving 
the gauge field A free. Using the new variables thus 
introduced, one can define a vector field n as follows. 
Let cr be the usual Pauli matrices. We then define 



i(lp*1p2 - IpllpZ) 

\^i\ 2 -\H 2 



where we demand that \^\ > everywhere, consider ^ 
normalised to unity and lim^^oo n = exists in order 
to obtain a map S 3 — * S 2 . Now the preimage of — 
forms a closed loop, the vortex core. 

The fact that A is left free, means that there is no non- 
trivial topology imposed on it and since in the vacuum of 
the Ginzburg-Landau model A is pure gauge, the mag- 
netic field energy can vanish in all cases. This in turn 
means that there is no longer a fourth-order derivative in 
the energy density (0]) and Derrick's theorem [l4[ states 
that no stable, topologically non-trivial solutions of the 
field equations with non-zero energy exist. Therefore, 
the collapse of the magnetic field must be somehow pre- 
vented in order to obtain stable topologically non-trivial 
configurations in the model. There are several physi- 
cal arguments that suggest there might exist physical 
processes that prevent the collapse, but here we follow 
the path set out by Ward, who used a geometrical argu- 
ment, by adding into the Lagrangian the term (we denote 
tft = tyj V 2 *)): 



AY 



(5) 



which makes the Ginzburg-Landau -Ward energy density 



£glw{x) 



-K 



i||D*|| 2 + I[ 

=£w 



V x A\\ 



(6) 



±K 2 \\^Dt>\\ 2 + V(iP u ip 2 ) 



and denoting for any subscript z: E z 
have the total energy 



J d 3 x£ z we finally 



Eglw = Ek + Ew + Eb + Ep. 



(7) 



The extra term, when the parameters k, rj — > oo, ensures 
that the model becomes exactly the Faddeev-Skyrmc 
model 



£fs = \\\dkn\\ 2 + \gFs\\n ■ d 3 n x <9 fc n|| 



(8) 



and therefore the model supports, at least asymp- 
totically, stable topologically non-trivial configurations; 
these solutions are called knot solitons due to their gen- 
eral shape. This limit of k — > oo was apparently first 
observed by Hindmarsh [ll[ , albeit in a slightly different 
context. 

Ward studied the question whether the solutions of the 
limiting model remain stable at finite values of re, rj. It 



was found, that if 77 = k 2 + 1, there are knot solitons 
already at k = 7.1. This is due to the fact that the extra 
term prevents the (total) collapse of the magnetic field, 
but only when k, rj are large enough: for smaller values, 
no solutions were found in Ref. [l 21 ] . although one was 
found in Ref. 

It was recently discovered independently by Babaev 
[l6| by using physical arguments and by Speight [lT], HI] 
by giving a rigorous mathematical proof, that the energy 
of the model has no topological lower bound, even when 
Cw is added, but instead for all values of the Hopf invari- 
ant, the energy can go to zero. Therefore any stable con- 
figurations found are necessarily only local minima of the 
energy; on the other hand, although the plain Ginzburg- 
Landau model does not seem to have any stable topolog- 
ically non-trivial configurations, they may only be very 
difficult to find due to very small attraction basin of said 
configurations, thus requiring very good initial guesses. 
Investigating such configurations in closely related mod- 
els may help in finding these initial configurations. 

In Ref. [12| Ward investigated only configurations, 
where 5 = 1 and rj = k 2 + 1. We will now present re- 
sults for the stable/unstable boundary rj{n) for g = 1 
and k € {6,20}. It was natural to start the investiga- 
tion from the values explored by Ward and expand the 
range of k in both directions; the limits of this range were 
eventually set by available computer capacity. 

It is also worth noting, that, as usual, Derrick's theo- 
rem provides a virial theorem for the model. Assume we 
have a solution of the field equations, "J, A and consider 
its energy density under uniform scaling of the coordi- 
nates x — > 7X: 

E G lw(i) = J d 3 x(£ K (jx) + £w{ix) + £b{jx) + £p{jx)) 

(9) 

which is the starting point of Derrick's theorem. Now, 
following the method used by Derrick, we get, after the 
change of integration variables x — > jx, 

E G lw{i) = lE K + j^Eb + jE w + j 3 E P (10) 
differentiating with respect to 7 we get 

-r 2 E B = E K + E W + 3j 2 E P . (11) 

The stability under scaling requires that this equation 
holds for 7=1, since otherwise some other size would 
be energetically more favorable. Thus we have a virial 
theorem: 



Ek 



3E P . 



(12) 



This must hold for all stable minimum energy configura- 
tions, regardless of the parameter values or value of Hopf 
invariant (topological charge) . For the remainder of the 
article, the energy is rescaled by 



E = 



Eglw 
4tt 2 V2 



3 



and we shall always use E for energy. This rescaling 
is motivated by Ref. and eases comparisons with 
Ref. [H. 

III. NUMERICAL METHODS 

We have discretized the system using single-step for- 
ward differences on a rectangular cubic lattice. Since 
the energy of the Ginzburg-Landau model is equal to 
that of the Abelian Higgs model, the discretization 
method is the standard (dropping the time-dependent 
part) used for lattice quantum field theories, as described 
in Refs. [li H3, EH . 

The use of single step in the finite differences approx- 
imation instead of some more sophisticated alternative 
with multiple points is simply a trade-off between speed 
and accuracy. The discretized equations are very long 
even with single step differences. This does not incur 
significant inaccuracy to the computation for two rea- 
sons. First, we are interested mainly in the existence of 
knotted solitons, which is not affected by the less accurate 
approximation - the exact values of the parameters where 
the transition from stable to unstable domain, however, 
do suffer from inaccuracies as shall be described later. 
Second, there is no accumulation of error during the it- 
erative process in the optimization algorithms employed. 

The term £ w was discretized in the same manner as 
the kinetic term: we denote fij = (<5j, <5|, Sj)v l , v l e 
{0, f } and find all gauge invariant discrete terms of the 
forms 

ip*(x + fi°)ij>i(x + (J,j) and 
^{x + n^ix + fi^), 

where v° + v 2 = v 1 + v 3 , J2i v < 4 and 

i\)\{x + [Xj)^i(x + [ij) and 

where X)j=o yl = There are 1 8 such terms, which are 
then multiplied by such constants that the sum of the 
multiplied terms has the correct continuum limit. This 
produces the most general single-step forward-differences 
discretization of £ w ■ 

Energy minimization was done using several differ- 
ent gradient-based optimization methods: steepest de- 
scents (SD), Fletcher-Reeves (FR) [22j and Polak-Ribiere 
(PR) [23| versions of the conjugate gradient method 
and Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi- 
Newton method 0, EE [H, [27j]. Of these a few simple 
speed tests were conducted, with the result that the FR 
method is usually fastest, but sometimes the SD method 
is faster due to the fact that the other methods spend 
too much time performing the line searches. (The effect 
of using less accurate line searches was not investigated.) 
All the methods are based on gradient directions and 



thus will only provide a local minimum. Some of the fi- 
nal configurations were subjected to simulated annealing 
in order to see how deep the minimum is. The annealing 
could not escape the minimum in a reasonable amount 
of time, making it reasonable to believe the minima are 
relatively deep. We did not have the computational re- 
sources to perform simulated annealing optimizations of 
all cases due to the extreme slowness of the algorithm. 

As a test of the accuracy and validity of the programs 
used, we reproduced the results of Ref. [l~2| section III. 
The results agree to within 5%, where our energies are al- 
ways higher; this confirms the correctness of the program 
and also gives some indication as to the accuracy com- 
pared to other methods. The accuracy of our method 
could be increased by using larger lattices and smaller 
lattice constants, but this kind of brute force approach 
would require excessive amounts of memory - we use more 
than half a terabyte at maximum - so some more sophis- 
ticated methods would be needed. 



IV. RESULTS 

The search for the boundary between stable and un- 
stable domains of (k^tj) was done as follows. First, an 
initial state was set up in a cubic 180 3 lattice with a lat- 
tice constant of 1/18. The initial configuration was con- 
structed so that none of the xy z-ax.es coincides with the 
axial symmetry of the soliton. This configuration was 
then minimized using one of the above algorithms and 
various combinations of (k, rf) to determine the rough 
shape of the boundary. This initial search is done in a 
small lattice in full knowledge that it may not be large 
enough to accurately distinguish between configurations 
which arc truly unstable and those that are unstabilized 
due to the large value of the lattice constant. Indeed, 
all the unstable systems at the boundary found in this 
initial search were later proven to be stable in a more ac- 
curate lattice. The search consisted of 72 computer runs, 
but due to the small lattice, used relatively little com- 
puter time and gave us the rough values around which 
to start searching for the boundary in a more accurate 
lattice. The number of these computationally much more 
expensive runs was 81. 

The values of (ac, rf) were then refined in lattices rang- 
ing from 360 3 to 600 3 . Some unstable configurations were 
also put into lattices of sizes up to 720 3 to support the 
conclusion that the instability is real and not caused by 
discretization effects. None of these were thus stabilised. 
The virial theorem Eq. (|12p was then checked for the sta- 
ble configurations at the boundary (i.e. for each n the 
stable configuration with lowest 77) to see if the config- 
uration really is a solution. Allowing for a 10% inaccu- 
racy, those that were within the tolerance were consid- 
ered solutions of the energy minimisation. These points 
are used to sketch the boundary of the stable configura- 
tions. Those that were outside the tolerance were further 
investigated. The reason for inaccuracy proved usually 
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to be due to the small physical dimensions of the compu- 
tational lattice: the knot suffers from pressure exerted by 
the edge of the box and cannot reach its preferred size. 
These configurations were therefore put into a larger lat- 
tice with the same lattice constant, minimized and the 
accuracy was checked again. This process was repeated 
as many times as necessary to achieve the desired ac- 
curacy - except for two cases as we will describe later. 
Whenever the accuracy was reached, the configuration 
was considered a solution and added to those used to 
sketch the boundary. For some cases the accuracy was 
simply a question of lattice constant; these were recom- 
puted with same physical dimensions but a smaller lattice 
constant in order to reach the desired accuracy. 

The exceptional cases where the process of putting into 
larger box until accuracy is achieved was not completed, 
were the pairs (k, 77) € {(10, 0.2), (20, 0.14)}. These are 
stable configurations, but the accuracy goal could not be 
achieved with the computational capability available due 
to the cubic growth of memory requirements of increasing 
the lattice size. 

After completing the above process, we select for each 
k the stable configurations with lowest 77; denote this 
value by 77™ 111 . The values (k, 77™ 111 ) are displayed as solid 
black circles in Fig. [T] together with a curve T){k) sketch- 
ing the approximate shape of the continuous boundary 
and yellow circles for largest unstable values of 77. Config- 
urations for different values of k are not always produced 
in a lattice of the same size, but despite that, the bound- 
ary curve fits rather well. All the stable dots in Fig. [T] 
are confirmed to be solutions by the virial theorem of 
Eq. p2[) . except the cases mentioned above: k <E {10, 20}. 
The boundary approaches y-axis as k — ► and x-axis as 
k — ► 00. The latter information is not very useful in 
constructing initial states for normal Ginzburg-Landau 
model, but the fact that the boundary seems to approach 
y-axis as well, might provide helpful insight and allow the 
construction of an initial state which could be used to find 
a topologically stable, non-trivial local energy minimum. 

Comparison of the energies of the final configurations 
reveals that the energy of the solution for (k, r) K ) is, 
within our numerical accuracy, independent of k. This 
is displayed in Fig. where the dots depict the ener- 
gies of the solutions, the heights of error bars are chosen 
according to how much the solution deviates from the 
virial theorem Eq. (fT2"]) and the solid horizontal line is 
the least-squares fit for the constant energy, neglecting 
the the anomalous cases k S {10,20}. Also, comparing 
this energy with the energies of the unstable configura- 
tions at the moment of loss of topology, shows that the 
unstable configurations always have lower energy than 
those at the boundary, giving even further support to our 
argument that the instability is real and not a numerical 
artefact. 

The search for the boundary also provides us, as a by- 
product, with information on the shape and size of the 
final configurations. All solutions have kept their initial 
orientation and overall shape, the only visible difference 
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6 8 10 12 14 16 18 20 

Figure 1: Boundary between stable and unstable 
regions: the solid black circles denote pairs of (k,i] k ), 
yellow circles denote the largest unstable values of 77 

and the curve, 77 = 45 7 - 65 ~ K -|- 0.5k~2 ; is a sketch of the 
boundary. 
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Figure 2: Energies (black discs) of the stable 
configurations closest to the boundary, the heights of 
the error bars show the inaccuracy of the solution as 
determined from the virial theorem and the horizontal 
line is the least-squares fit of a constant energy. 



between the final configurations is the apparent decrease 
of the size of the torus with growing k, as shown in Fig. [3] 
for the cases k€{6,8,12}. It remains an open question 
whether the final toroidal configuration obtained from 
this initial configuration would shrink to zero as n — > 00 
because we were unable to follow the boundary above 
k = 12. 

Looking at the four terms of the energy function of 
these final configurations reveals more details of the in- 
terplay between the various terms. As can be seen from 
Table HI Ek increases with increasing k, but Ew de- 
creases. This is expected since on the limit k — > 00, we 
must have Ew — > 0. Since it was found in [l3[ that the 
magnetic field always approaches to zero for unstable sys- 
tems in non-modified Ginzburg-Landau model, it is in- 
teresting to note that there seems to be no trace of this 
here: the magnetic energy does not change appreciably. 
The variation in the potential energy is also negligibly 
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Figure 3: Initial configuration (top left) and three 
resulting final configurations: n = 6 (top right), k = 8 
(bottom left) and k = 12 (bottom right). The coloured 
region is the equator of n (i.e. isosurface where n 3 = 0) 
and the coloring corresponds to the longitude of n. 
(Colors available on-line.) 

small. 



K 


E K 


E B 


Ew 


E P 


6 


23.740 


32.111 


7.6571 


0.46053 


8 


26.739 


31.808 


5.9745 


0.46335 


12 


31.580 


32.334 


2.5649 


0.70555 



Table I: Terms of the energy function of the final 
configurations of k € {6, 8, 12}. 



each k below which the minimum vanishes. It remains 
open whether there still is a minimum but our initial con- 
figuration has simply moved "closer" to the global min- 
imum of zero so that the gradient-based algorithms can 
no longer reach it. Also, there can be other local minima. 
To explore these possibilities further it is required to use 
either a set of very different initial configurations or an 
algorithm which can explore a wide region of the config- 
uration space starting from a single initial configuration 
such as the genetic algorithm (which has the downside of 
being able to escape the local minima and thus ending 
up in the trivial global minimum). 

Strikingly, the energy was shown to be constant along 
the boundary. This information, along with the fact 
that higher values of r] are required for lower values of 
k, might provide a way to construct an initial configu- 
ration in the ordinary two-component Ginzburg-Landau 
model, which, under minimisation of energy, would lead 
to a non-zero local minimum. The procedure would, how- 
ever, require further insight into how the various terms 
of the energy functional behave when k decreases; our 
numerical scheme was not designed for this and as such, 
appears not to be accurate enough to provide this infor- 
mation. In contrast to its energy, the size of the minimum 
energy configuration decreases as k increases. This re- 
quires progressively smaller values of the lattice constant 
and thus significantly different numerical approach than 
the simple, but very large (recall that we used over half a 
terabyte) lattices used here. Still, the possibility remains 
of further research in both smaller and larger values of n, 
but as they are not addressable by the framework used 
here, it falls outside the scope of this paper. 
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V. CONCLUSIONS 

We have studied the existence of local minima in a 
modified two-component Ginzburg-Landau model and 
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where the local minima become unstable. It was found 
that local minima exist for a wide range of values of the 
parameter k, but that there is a limiting value of r\ for 
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